function   results =  solve_v(  model_sol,  par  )
% This file sovles the bond price v for both risky bonds and the safe bonds. 
% Use false time derivative method. 
% alg_control is a subfield of par.  It contains information about whether to plot and show the progress.

options = optimset('TolX',1e-30);
lambda_star = fzero( par.mu_lambda_fun,  (par.lambda_H+par.lambda_L)/2 , options );
par.lambda_star = lambda_star;
tau = par.tau;   % bond maturity for credit spread pricing

% === Think about this file as a function. Parameters include the raw model S, the maximum of maturity, the current maturity, and loss given default.  
N_w = par.N_w;  
N_lambda = par.N_lambda; 
w_grid = par.w_grid;  lambda_grid=par.lambda_grid;
w_matrix = model_sol.w_matrix; 
lambda_matrix = model_sol.lambda_matrix; 

alg_control = par.alg_control;
if(  ~ (isfield( par.alg_control, 'plot' ) && par.alg_control.plot) )
    alg_control.plot = false;
end

%% Step 0:  Find the Boudary Conditions at w=0 and w=1
% ***** Part 0: find the critical boundary lambda*, and the risk-free rates at the boundary w=0 and w=1  ***** 
w_dim = kron(ones(N_lambda,1), w_grid(2:(N_w-1)) );
lambda_dim = kron(lambda_grid, ones(N_w-2,1));
fit_a_function = @(values_matrix)  scatteredInterpolant(w_dim,    lambda_dim ,  reshape( values_matrix(2:(N_w-1),:) , (N_w-2)*N_lambda, 1 ),   'linear',   'linear')  ;   % note: linear extrapolation at the boundary

% ***** Part 1: boundary conditions at w=0 and w=1***** 
% Note: we have a jump component in the equation. Therefore, we need to solve for jump problems. 
w=0;   
options.safe_bonds =  false;     options.plot=false;
result = bond_pricing_solve_boundary_of_w( w, par,  model_sol, options );
lambda_vec = result.lambda_vec;    v0_vec = result.v_vec;    v0_fun_lambda = result.v_fun_lambda;

options.safe_bonds = true; 
result = bond_pricing_solve_boundary_of_w( w, par,  model_sol, options );
v0_safe_bond_vec = result.v_fun_lambda(lambda_vec) ;

% risky bond pricing at w=1
w=1;   
options.safe_bonds = false;   
result = bond_pricing_solve_boundary_of_w( w, par,  model_sol, options );
v1_vec = result.v_fun_lambda(lambda_vec) ;    v1_fun_lambda = result.v_fun_lambda;

options.safe_bonds = true; 
result = bond_pricing_solve_boundary_of_w( w, par,  model_sol, options );
v1_safe_bond_vec = result.v_fun_lambda(lambda_vec) ;

% comparison of solutions
if(alg_control.plot)
    figure; 
    subplot(1,2,1); plot( lambda_vec, v0_vec,  lambda_vec, v0_safe_bond_vec, '--', 'LineWidth',2 ); xline(lambda_star, 'k:'); xlabel('\lambda'); ylabel('bond price at w=0');  legend('risky bonds', 'safe bonds');
    subplot(1,2,2); plot( lambda_vec, v1_vec,  lambda_vec, v1_safe_bond_vec, '--', 'LineWidth',2  ); xline(lambda_star, 'r-.'); xlabel('\lambda'); ylabel('bond price at w=1');   legend('risky bonds', 'safe bonds');
    pause(2); close all;
end

% % ***************************************************************************************
% %                   Method 1: Solve from False Time Iterations 
%% **************************************************************************************************
h = 0.05;
% Use the boundary solutions to initialize the problem. 
w_init = [  zeros(size(v0_vec));     ones(size(v1_vec))    ];    %;    w_grid
lambda_init = [  lambda_vec;  lambda_vec     ];   %;   lambda_star*ones(size(w_grid))
v_init =  scatteredInterpolant(w_init,lambda_init,  [ v0_vec; v1_vec ]   ,'linear','linear');   %;  v_lambda_star

for( round = 1:2)   % round=1 solves the safe bonds, while round=2 solves the risky bonds
    if(round==1)
        hat_kappa0 = 0;   % safe bonds
    else
        hat_kappa0 = par.hat_kappa0;  
    end

    w_dim = kron(ones(N_lambda,1), w_grid );
    lambda_dim = kron(lambda_grid, ones(N_w,1));
    fit_a_function_full = @(values_matrix)  scatteredInterpolant(w_dim,    lambda_dim ,  reshape( values_matrix , N_w*N_lambda, 1 ),   'linear',   'linear')  ;   % note: linear extrapolation at the boundary
    v_new = v_init; 
    v_last = @(w, lambda) v_init(w, lambda) + 1 ;
    %  plot(   w_grid,  v_init(w_matrix,lambda_matrix)    )
    %  plot(   lambda_grid,  v_init(w_matrix,lambda_matrix)    )
    iter_count = 0; max_it = 500;      diff_threshold = 0.02;   diff=1;
    if((par.model_name=="benchmark"))
        diff_threshold = 0.001;
    end
    if( par.model_name=="behavioral" )
        diff_threshold = 0.005;
    end
    diff_vec = [];
    while( iter_count < max_it  &&  diff > diff_threshold )
        iter_count = iter_count+1; 
        diff_matrix= v_new(w_dim, lambda_dim) - v_last(w_dim, lambda_dim);
        diff = sum(sum( sum(  abs(  diff_matrix   )  ) ));
        if(  isfield( alg_control, 'show_progress' ) && alg_control.show_progress && round==2)
            if( mod( iter_count, 20 )==0 )
                disp( ['Solving credit spread: Iteration ', num2str(iter_count), ',  diff=', num2str(diff) ] )
            end
        end
        diff_vec = [diff_vec; diff];
        v_last = v_new; 
        v_last_matrix = v_last(w_matrix, lambda_matrix);
        % method: Runge Kutta 4. Later discovered that improving integration accuracy does not help much
        k1 = h*time_derivative_v(  v_last, model_sol,  par, hat_kappa0 );
    %     k2 = h*time_derivative_v(  fit_a_function( v_last_matrix - k1/2 )  , model_sol,  par, hat_kappa0);
    %     k3 = h*time_derivative_v(  fit_a_function(  v_last_matrix - k2/2 ) , model_sol,  par, hat_kappa0);
    %     k4 = h*time_derivative_v(  fit_a_function(  v_last_matrix - k3 )  ,   model_sol,     par, hat_kappa0);
    %     v_new_matrix = v_last_matrix - 1/6*( k1 + 2*k2+2*k3 + k4 );
        v_new_matrix = v_last_matrix - k1;
        % Generate the fitted functions
%         v_new = fit_a_function_full( v_new_matrix );
%         [~, v_new_matrix] = time_derivative_v(  v_new, model_sol,  par, hat_kappa0 ) ;
        v_new_matrix(1,:) =  v0_fun_lambda(par.lambda_grid) ;   % boundary condition at w=0
        v_new_matrix(end,:) = v1_fun_lambda(par.lambda_grid) ;   % boundary condition at w=1   
        
        % Beginning of smoothing -- Sometimes the smoothing is useful.
        v = v_new_matrix;
        for( j = 1:par.N_lambda)
            v(:,j) = smooth( w_grid, v(:,j) );
        end
        for( i = 1:par.N_w)
            v(i,:) = smooth( lambda_grid, v(i,:));
        end
%         if(round==2)
%             v =  min( v, v_safe(w_matrix, lambda_matrix)  ) ;   % note: risky bonds are less expensive than risk-free bonds.
%         end
        % End of smoothing
        v_new = fit_a_function( v );
        
        if( mod( iter_count, 1 )==0 && alg_control.plot )
            figure(1);
            subplot(1,2,1);  plot(  w_grid, v_new_matrix(:,5) ); xlabel('w'); ylabel('v');
            subplot(1,2,2);  plot(  lambda_grid, v_new_matrix(20,:) );  xlabel('lambda');  ylabel('v');
            pause(0.2); 
        end
        if( iter_count>=2 && diff_vec(end)>diff_vec(end-1) )
            h = h/2;
        end
    end
    close all;

    % Plot the solution
    if(alg_control.plot)
        figure;  
        subplot(1,2,1); plot( par.w_grid, v_new_matrix );   xlabel('banker wealth w');  ylabel('bond price'); 
        subplot(1,2,2); plot( par.lambda_grid, v_new_matrix );   xlabel('lambda');  ylabel('bond price'); 
        pause(2);
    end

    if(round==1)
        v_safe = v_new;
    else
        v_risky = v_new;
    end

    v_init = v_new;
end

credit_spread =  ( log( v_safe(w_matrix, lambda_matrix)  ) -  log(  v_risky(w_matrix, lambda_matrix)  )  ) / tau;

if(alg_control.plot)
    figure;  
    subplot(1,2,1); plot( par.w_grid, credit_spread );   xlabel('banker wealth w');  ylabel('credit spread'); 
    subplot(1,2,2); plot( par.lambda_grid, credit_spread );   xlabel('lambda');  ylabel('credit spread'); 
    pause(0.5);
end

results.v_safe = v_safe;
results.v_risky = v_risky;
results.credit_spread = credit_spread;






